#define NR 6
#define NC 4
#include <stdio.h>
#include "nrutil.h"
#include "nr.h"

main()
{

    static double c[6][4]={
			  {0.,0.,0.,0.},
			  {0.,0.,0.,0.},
			  {1.,1.,1.,1.},
			  {-25.,-25.,25.,25.},
			  {15.,-15.,-15.,15.},
			  {0.,0.,0.,0.}
			};

    double wmax,wmin,**a,**u,*w,**v,*b,*x;
    int i,j;



    a=convert_matrix(&c[0][0],1,NR,1,NC);
    u=matrix(1,NR,1,NC);
    v=matrix(1,NC,1,NC);
    w=vector(1,NC);
    b=vector(1,NR);
    x=vector(1,NC);
    for(i=1;i<=NR;i++) {
      for (j=1;j<=NC;j++) {
	u[i][j]=a[i][j];
      }
    }

    svdcmp(u,NR,NC,w,v);

    printf("singular value=");
    for(i=1; i<=NC; i++)	printf(" %f ",w[i]);
    printf("\n");

    for(i=1; i<=NC; i++)
	printf("u[%d]= %f %f %f %f %f %f\n",i,
			u[1][i],u[2][i],u[3][i],u[4][i],u[5][i],u[6][i]);

    for(i=1; i<=NC; i++)
	printf("v[%d]= %f %f %f %f\n",i,
	       v[i][1],v[i][2],v[i][3],v[i][4]);

wmax=0.0;
for(j=1;j<=NC;j++) if (w[j] > wmax) wmax=w[j];
wmin=wmax*1.0e-6;
for(j=1;j<=NC;j++) if (w[j] < wmin) w[j]=0.0;

for(j=1;j<=NR;j++) b[j]=(double)j;

svbksb(u,w,v,NR,NC,b,x);

printf("x=   ");
for(i=1;i<=NC;i++) printf(" %f",x[i]);
printf("\n");

}
